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Motivated by experiments on Josephson junction arrays in a magnetic field and ultracold inter- 
acting atoms in an optical lattice in the presence of a 'synthetic' orbital magnetic fields, we study 
the "fully frustrated" Bose-Hubbard model and quantum XY model with half a flux quantum per 
lattice plaquette. Using Monte Carlo simulations and the density matrix renormalization group 
method, we show that these kinetically frustrated boson models admit three phases at integer fill- 
ing: a weakly interacting chiral superfluid phase with staggered loop currents which spontaneously 
break time-reversal symmetry, a conventional Mott insulator at strong coupling, and a remarkable 
"chiral Mott insulator" (CMI) with staggered loop currents sandwiched between them at interme- 
diate correlation. We discuss how the CMI state may be viewed as an exciton condensate or a 
vortex supersolid, study a Jastrow variational wavefunction which captures its correlations, present 
results for the boson momentum distribution across the phase diagram, and consider various exper- 
imental implications of our phase diagram. Finally, we consider generalizations to a staggered flux 
Bose-Hubbard model and a two-dimensional (2D) version of the CMI in weakly coupled ladders. 



I. INTRODUCTION 

The effect of frustration in generating unusual states 
of matter such as fractional quantum Hall fluids or quan- 
tum spin liquids is an important and recurring theme in 
the physics of condensed matter systems Recently, re- 
search in the field of ultracold atomic gases has begun to 
explore this area, spurred on by the creation of artificial 
gauge fields using Raman transitions in systems of cold 
atoms These gauge fields can be used to thread fluxes 
through the plaquettes of optical lattices giving rise to 
"kinetic frustration" by producing multiple minima in 
the band dispersion and frustrating simple Bose conden- 
sation into a single non-degenerate minimum. Similarly, 
time-dependent shaking of the optical lattice^ or pop- 
ulating higher bands of an optical lattice^ can be used 
to control the sign of the hopping amplitude in an op- 
tical lattice, again leading to such "kinetic frustration" . 
For bosonic atoms with weak repulsion, such kinetic frus- 
tration gets resolved in a manner such that the resulting 
superfluid state can have a broken symmetry correspond- 
ing to picking out a particular linear combination of the 
different minimai^r— Increasing the strength of the in- 
teractions at commensurate filling can be expected to 
eventually yield a Mott insulator (MI) with the motion 
of the bosons quenched, which thus renders the kinetic 
frustration ineffective. In a synthetic flux and at strong 
coupling, the fully gapped MI is identical to the one ex- 
pected for the same lattice without a frustrating flux per 
plaquette^ this simply means that at strong coupling, 
we can adiabatically remove the flux without encounter- 
ing a quantum phase transition. However, there could 



exist a state intermediate to the superfluid and the MI 
described above for which charge motion has been sup- 
pressed enough to open up a gap but not to restore the 
broken symmetry associated with frustration. Such a 
'weak' Mott insulating state, while it is globally incom- 
pressible, supports significant local boson number fluc- 
tuations, so bosons can 'sense' the magnetic flux on a 
plaquette. In a recent paper, we have found numerical 
evidence for the existence of such a remarkable interme- 
diate state in frustrated two-leg ladders of bosons for the 
so-called Fully Frustrated Bose Hubbard (FFBH) model 
which has half-a-flux-quantum per plaquette. We call 
this state a "chiral Mott insulator" (CMI) since it is fully 
gapped due to boson-boson interactions, exactly like an 
ordinary Mott insulator, and in addition possesses chi- 
ral order associated with the spontaneously broken time- 
reversal symmetry arising from resolving the kinetic frus- 
tration. The superfluid state of this system also possesses 
this chiral order and we thus dub it a chiral superfluid 
(CSF)J^ Other recent studies have also focussed on a 
variety of such exotic states driven by "ring-exchange" 
interactions,— ~— which again arise in a similar fashion 
driven by virtual charge fluctuations in a Mott insulator. 

In this paper, we discuss further details of our work on 
this FFBH ladder model, and its close cousin, the fully 
frustrated quantum XY model to which it reduces at high 
filling factors. We also discuss how one might stabilize 
such a Mott insulator in higher dimensions by considering 
weakly coupled ladders. Such a CMI may also be viewed 
as a bosonic Mott insulating version of the staggered loop 
current states^— studied in the context of pseudogap 
physics in the high temperature cuprate superconductors. 
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Classical analogues of the CMI and CSF states have 
been studied in the paster— The simplest classical model 
displaying analogous phases is the fully frustrated XY 
model in two dimensions! 23 ' 24 At small but finite temper- 
ature, this model has a phase with algebraic U(l) order 
for the spins along with a staggered pattern of vorticity 
associated with each plaquette corresponding to broken 
Z2 symmetry. This is the analogue of the CSF phase. 
As the temperature is increased, the U(l) symmetry is 
restored while the Z 2 symmetry continues to be broken 
in a state that is the analogue of the CMI. Upon fur- 
ther increasing the temperature, the Zi order is restored 
yielding a completely disordered state, which is the ana- 
logue of the featureless CMI. Analogs of these phases 
have also been found in another classical models which 
will be discussed late r 43 ' 44 . 

The CSF phase has been studied in a variety of frus- 
trated quantum models of bosons but without the strong 
correlations required to obtain the CMI phase.— ~— The 
effect of correlations in conjunction with frustration has 
been investigated for models of fcrmions and quantum 
spinsi 20 ' 22 ' 29 " — In both cases, analogs of the CSF state 
are obtained as staggered current (for fermions) or gap- 
less chiral (for spins) states. An analog of the CMI state 
has been found for fermions in a ladder models For 
spins, the analog of the CMI would be a spin gapped 
state with vector chiral order, which has been proposed 
for easy-plane frustrated magnets.— Studies on micro- 
scopic spin models for a CMI-like state suggest that it 
coexists with either dimer order (for half integer spin) or 
topological order (for integer spin) . 38 ' 39 

In this paper, we study the Bose-Hubbard ladder of 
ref. [ID using a variety of different techniques to elucidate 
the nature and properties of the CSF and CMI. We em- 
ploy mean field theory, mapping onto an effective clas- 
sical model and a variational Monte-Carlo method for 
this purpose. We also provide physical pictures for the 
CMI phase as a supersolid of vortices and a condensate of 
neutral excitons. The outline of the paper is as follows: 
In section II, we introduce the microscopic model with a 
summary of results obtained from numerics in our previ- 
ous work. In section III, we perform a mean-field calcula- 
tion which provides a good description of CSF phase but 
is unable to describe the CMI phase. Then in section IV, 
we write down a rotor model for our Hamiltonian which 
allows us to map it onto a model, which can be studied us- 
ing classical Monte-Carlo simulations. Numerics on this 
model have been performed in ref. [TH and show the exis- 
tence of the CMI. This resulting model can be used for 
comparison to the classical frustrated models mentioned 
earlier. In section V, we provide details of our density 
matrix renormalization group (DMRG) studies on the 
Bose-Hubbard model, which also reveals a CMI phase. In 
section VI, we provide physical pictures of the CMI as a 
supersolid of vortices or a condensate of neutral excitons 
while in section VII, we perform a variational Monte- 
Carlo study of a candidate Jastrow wavefunction for the 
CMI state, showing that it correctly captures the essen- 
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FIG. 1. Sketch of the fully frustrated Bose Hubbard two-leg 
ladder. The opposite signs of the hopping on the a and b 
legs correspond to a 7r-flux threading each square plaquette 
and frustrating the boson kinetic energy. Also depicted at the 
spontaneously generated staggered loop currents generated in 
the superfluid and chiral Mott insulator phases of this ladder 
model. 



tial correlations of the CMI state. Section VIII provides 
a short description of experiments that can be performed 
to detect the CMI state in Josephson junction arrays or 
cold atom systems. Section IX discusses generalizations 
to a staggered flux model and to higher dimensions, and 
we conclude in Section X with a summary. 



II. MODEL AND PHASE DIAGRAM 

The Hamiltonian for the frustrated ladder as shown in 
Fig. [1] can be written as 

X X 
X X 

where a x and b x are bosonic operators on each of the two 
legs of the ladder whose sites are labelled by x. n a ^ x and 
Ub tX ar e the corresponding occupation numbers. tj_ is the 
hopping amplitude between the legs and U is the onsitc 
two-body interaction strength. 

In our previous work we studied the ground state of the 
Eqn. [1] at integer filling using two independent numeri- 
cal methods such as classical Monte Carlo techniques^ 3 , 
and density matrix renormalization group(DMRG). The 
complete phase diagram of this model is shown in Fig(|5|) 
and Fig(jn]) and described in detail in the corresponding 
sections. The important results obtained by the above 
two analysis are qualitatively similar. As a result of the 
competition between the onsite interaction ([/), intra- 
chain hopping( t) and interchain hopping (t±), we ob- 
tain three different quantum phases: the CSF, the CMI 
and a regular MI . When the Hubbard repulsion (U) is 
small the system exhibits a gaplcss SF phase with a fi- 
nite loop current order in each plaquette. For intermedi- 
ate values of U the system undergoes a transition from 
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CSF to CMI phase which possesses finite charge gap and 
also exhibits staggered loop current and spontaneously 
breaks time reversal symmetry. A further increase of the 
value of U breaks the loop current order in the system 
and the system undergoes a transition to the regular MI 
phase. Using different scaling properties of the observ- 
ables obtained by the numerical simulations we found 
that the SF to CMI phase transition is of Bcrczinskii- 
Kosterlitz-Thouless (BKT)iS universality class and the 
transition from CMI to MI is of Ising universality class. 
Apart from the numerical analysis of this model we anal- 
yse the system using several analytical methods. We ex- 
plain the numerical and analytical methods used to un- 
derstand the ground state properties of the Hamiltonian 
given in Eqn. [TJ and describe the properties of the differ- 
ent phases, in the following sections. 



III. APPROXIMATE ANALYSES OF THE FFBH 
LADDER MODEL 

A. Single particle condensate wavefunction 

We begin by motivating the form of the single particle 
condensate wavefunction for our system in the weakly 
interacting limit. To do this we first consider the sin- 
gle particle dispersion obtained from Eqn. [TJ i.e. setting 
[7 = 0. If we also set t± = 0, the dispersion is as shown 
in the upper panel of Fig. [2] with two bands inverted with 
respect to each other and intersecting at k = ±7r. The 
band with a minimum at k — corresponds to the a leg 
while the one with a maximum at k = corresponds to 
the b leg. With t± ^ 0, the degeneracies at the points of 
intersection are lifted resulting in two bands with a gap 
as shown in the lower panel of Fig. [2j The two minima at 
k = and k = n originate from the bands corresponding 
to each of the two legs with the k = minimum cor- 
responding to bosons localized in the a leg and k = it 
minimum to bosons in the b leg. A single particle con- 
densate wavefunction for U — can thus be any linear 
superposition of the states at the two minima. 

For U > at a filling of one boson per site, the wave- 
function has the form 

|*> = ^c <9 (|* = 0) + e < *|* = 7r)) 1 (2) 

where \k = 0) and \k = ir) correspond to the states at 
the minima k = and k — tt respectively and 9 and <fi 
are phases. This form of comes about because the 
onsite repulsion discourages double occupancy so that 
with one boson per site, the system favors distributing 
the bosons equally over both legs. Thus, the single parti- 
cle condensate wavefucntion has equal amplitude in the 
states \k = 0) and \k = ir) with relative phase c\> and 
global phase 9. Minimizing the energy of states of the 
form given in Eqn. [5]yields <fi = ±7r/2. Thus, the form of 




FIG. 2. (Top) The single particle dispersion obtained from 
Eqn. Q] with tj_ = 0. The two legs are now completely de- 
coupled from each other. For t > 0, the lower (upper) band 
corresponds to the bosons being in the a (6) leg as is shown 
by the solid (dashed) line.. The two bands are degenerate at 
k = ±7r/2. (Bottom) For t± ^ 0, the degeneracies between 
the bands are lifted resulting in a gap between the two bands. 
There are now two minima in the lower band at k = and 
7T. The k = (n) minimum originates from the lower (upper) 
band of the two decoupled chains and thus corresponds to the 
particles being localized mostly in the a (b) legs. 

the single particle wavefunction is 

^) = ±e^(\k = 0) + e ±i ^ 2 \k = 7r)), (3) 

which has a global U(l) symmetry corresponding to 9, 
which acts as the supernuid parameter and a Z 2 symme- 
try due to the ± sign in the relative phase corresponding 
to a particular pattern of staggered current loops (chi- 
rality). In 1 + 1 D, like for our system, there is no long 
range U(l) order but at best algebraically decaying cor- 
relations corresponding to quasi-long range order of the 
supernuid. However, there can exist true long range Ising 
order. The state with coexisting algebraic supernuid or- 
der and long range chiral order is the chiral supernuid. 
The chiral Mott insulator corresponds to short range su- 
pernuid correlations (accompanied by a charge gap) ex- 
isting along with long range chiral order whereas in the 
Mott insulator both supernuid and chiral order are com- 
pletely absent. 

B. Mean field theory of the chiral superfluid 

In this section we refine the above discussion and 
present the mean-field theory description of the chiral 
superfluid phase in ([1}. Ignoring the interaction term 
U , we obtain the single particle dispersion ±£7fe with 
Ek = \/tj~+ (2t cos fc) 2 and introducing a chemical po- 
tential fi, the two bands have dispersion ±Ek — \i. Of 
these, the lower band (with dispersion — Ek — /x) has 
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two minima at k = 0, it. The eigenvector for this band 
(uk,v k ) T is 



1 



v/1 



9l 



2t cos k 
E k 



Vk = 



9k 



u k v k 



t± 

2E k 



1 2tcosk^ 

2 E k 



(4) 

(5) 
(6) 



where g k = (Ek — 2t) / 1± . Note that we can set u k +Tv = v k 
and Vfc+Tr = u k which means, in particular, that u T = vq 
and = uq. 

We can rotate to a new basis, where at (ft].) creates 
quasiparticles in the lower (higher) band at momentum 
k, via 



k - 



u k -v k 
v k u k 



ft,, 



(7) 



in which the single particle Hamiltonian is diagonal and 
given by 



]T [(-E k - n)ala k + (E k - ^fifa) 

k 



(8) 



same for both, which we call ip. 
global phase rotation, 



This gives us up to a 



(a x ) = i)[u ± iv (-l) x ] 
(b x )=i>[v Q ±iu Q (-lf] 



(11) 
(12) 



This leads to a spatially uniform boson density and bond 
energy, but to a spatially varying bond current density 
which forms a staggered pattern, much like a vortex- 
antivortex crystal. Specifically, the bond currents are 
given by 



3x,x+l — tt{ a x a x+l 



l x a x+ll 



Jx,x-\-l 



Jx = 



T4# 2 Uo « (-ir 

+ it ( b l b x+i - b lK+i) 

±4ti/j 2 u v (-l) x 

-it±{a>th - b l a x) 
T2t x i?(ul-vl)(-iy. 



(13) 

(14) 
(15) 



Note that current conservation at each lattice point works 
out if we use the fact, coming from the single particle dis- 
persion analysis, that 4tuoVo = t±(ug — Vq). The two pos- 
sible choices for the sign of each bond current correspond 
to two distinct current order patterns which are related 
to one another by time-reversal or by a unit translation. 



Reintroducing the local repulsive interaction iJi n t = 
\(r? a + n%) at each site, which is responsible for eventu- 
ally driving the 2-leg ladder into a Mott insulator state, 
we get back the Hamiltonian in Eqn. [T] We focus on the 
low energy physics and thus ignore any effects of the up- 
per band. Further, at the mean-field level we can focus 
only on the lowest energy k = 0, 7r modes in the low en- 
ergy band, and write the Hamiltonian in terms of these 
modes. This leads to 



HZ> = (-E - /x) 2 «J< 



-U(u 4 



i = 0,7T 

,4 



at at a, a. 



1 — 0.7T 



+Wulvla\ala v a 

+2Uulvl(a\ala„a^ + a\ala a a Q ) (9) 

If we were to condense the 04 bosons with (ctj) = tpi 
(where i = 0,7r), we get the mean held energy 



E 



i=0, 7T 



+8Uu 2 v%\ip \ 2 \<p n \ 



J2 \<p* 

1— 0,7T 



u(4 + v i ) 



where the hnal term describes "Umklapp effects" which 
transfers a pair of bosons from one minimum to the other. 
Minimizing the interaction energy with respect to the 
phase difference between tpo and (p^ gives a value of ±n/2 
phase for this quantity. Minimizing with respect to the 
amplitude of the two condensates we find it to be the 



C. Mean field theory of the Mott transition 

We next turn to the Mott insulating state induced by 
strong repulsion in the CSF state. We describe here the 
mean field theory of this Mott transition, which we then 
argue to be inadequate to describe the physics. To ob- 
tain the mean field CSF to insulator phase boundary, we 
consider the single site Hamiltonian (motivated by our 
previous discussion of the CSF), 



h: 



Slllfi 

f 



— 2ti/j((u — ivo)a' + (uo + iv )a ) 

t±ip((vo + iu )a^ + (v - iu )a ) 
U 



— a'a'a a 
2 



(16) 



Setting t and t± to we obtain a single site Hamiltonian, 
whose ground state is a true number state |n) given by 
U(n — 1) < yU < Un. Turning on a small nonzero ip leads 
to the first order corrected wavefunction 



where 



ip(r - is)y/n + 1 . , ip(r + is)y/ri . 

-|n+l)H — — In— 1) 



Un — fi 



fi - U(n - 1) 



r = (2tu + txv ), 
s = (2tvo — t±uo). 



(17) 

(18) 
(19) 



The resulting linearized self-consistent equation for ip 
leads to the phase boundary 



1 



1 



y/At 2 + tl fi - U(n - 1) Un - fi 



(20) 
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with U(n — 1) < fx < Un. 

This result for the SF-MI phase boundary is exactly 
the same as in the usual mea n field SF -MI phase diagram 
provided we replace zt — > \/4t 2 + t\ in the conventional 
result (where z is the coordination number), so that the 
Mott transition happens when 

r r7r— flux 

c.ladder \ 

where « 5.83 at a filling of one boson per site. 

We end by noting that this mean field theory does not 
have explicit current density wave and superfluidity as 
two independent order parameters, since the current is 
simply a product of the Bose condensate at two differ- 
ent sites. Thus the superfluidity and any time-reversal 
breaking orders vanish together at a single SF-MI transi- 
tion. However, such a vanishing of two completely differ- 
ent orders at a single continuous Mott transition is not 
expected to be generic and is Landau-forbidden. This 
encourages us to pursue careful numerical studies of this 
model which enable us to go beyond this simple mean 



field theory. 

IV. CLASSICAL XY MODEL 

As noted above, mean field theory fails to give us 
the complete picture of the FFBH model. In order to 
go beyond this approach, we map the quantum FFXY 
model (which is a good description of the FFBH model 
at large integer filling) to an effective classical model on a 
space-time lattice by stacking of ladders atop each other 
in the imaginary space time. To do this we construct 
a rotor model from the fully frustrated Bose Hubbard 
model in Eqn. [T] - such a rotor model ignores amplitude 
fluctuations and is expected to be a valid effective de- 
scription at large fillings. Setting a x ~ exp(— itp*) and 
b x ~ exp(— iip x ) and replacing the number operators n a ' b 
by angular momentum operators C a,b which cause fluctu- 
ations in the phases (angles) (p a ' b we obtain the partition 
function in terms of the classical action in one higher 
dimension as 

Z = E e" s i +1 M (22) 

where 



5"cl +1 = - £ [ J H C0S (^+l,r - <Px,r) - J\\ COs{(p b x+1>T - tp XiT ) + J± COsfo£ lT - <4,t)] 



Jr E [ C0S «r+l - <P%,r) + C0S «r+l ~ ¥>x,r)] > 



(23) 



with 2et = Ju, 2et± = J±, and I /ell = J T (see Appendix 
A for a detailed derivation). We see that this has the 
form of a XY model. 



H 



XY 



2J Js COS (tpi - tfi+s) . 



(24) 



Here the classical variable ipi corresponds to the boson 
phases and (i,i+S) are the nearest neighbours along the 
space-time direction 6. The couplings Jg take the values 
±J|I on the two legs, Jj_ on the rungs linking the two 
layers, and J T in the imaginary time direction. In order 
to get the properties of the quantum model at a fixed in- 
verse temperature /3t, we must take the 'time'-continuum 
limit of e — > 0, sending Jy — > 0, J± — > 0, and J T — > oo, 
while keeping fixed J±/J\\ = tj_/i and J T J\\ = 2t/U. The 
inverse temperature jit is then given by etL T and thus 
depends on the chosen value of et (which must be taken 
to be very small) and the size of the simulation cell in 
the 'time'-direction. We set e = I / \JlUi which leads to 

J,| =J T = yWU- 

This effective square lattice bilayer XY Hamiltonian 
for the FFBH model can be simulated using a classical 



Monte Carlo algorithm. The Monte Carlo simulation was 
performed using the Metropolis algorithm with selective 
energy conserving moves to increase the acceptance rates 
at large coupling. About I0 6 — 10 7 steps were used to 
calculate the equilibrium values of the various quantities 
at the largest system sizes. 

The ordinary Mott insulator of the FFBH model cor- 
responds to a fully disordered paramagnetic state of the 
effective bilayer classical XY model. Starting from the 
supcrfluid phase at small 1/ J T , we detect the vanishing 
of superfluidity in this model by calculating the helic- 
ity modulus T, which corresponds to a response to an 
infinitesimal phase twist in the spatial direction. Explic- 

itly, 



r= l J 2F 



2 



*->o 



where the free energy 



F = -log ]T e~ s - , 

{Vx, T } 

and $ is the flux twist along the || direction. 



(25) 



(26) 




FIG. 3. Phase diagram of the effective classical model ob- 
tained in ref. from classical Monte Carlo simulations of 
Eqn. 1241 which is shown in with = J T , obtained from vari- 
ational Monte Carlo simulations. 



As seen from Fig. 2] T appears to drop precipitously 
in the range 1/J T ~ 0.8 — 1.1, with the size dependence 
of r indicative of a jump in the thermodynamic limit, as 
expected at a BKT transition. However, for finite-sized 
systems, this jump is rounded off and one has to use 
finite-size scaling to precisely locate the transition. The 
relevant equation which governs the behavior of the he- 
licity modulus at the transition point is obtained from in- 
tegrating the Kosterlitz-Thouless renormalization group 
equations, and yield o 24 ' 41 

™- A ( 1+ ki*kc)- (27) 

where the constant A has the universal value of 2/tt for 
the BKT transition while C is a non-universal constant 
and the classical system is of size L x L x 2 (with the 
factor of 2 for our bilayer). T(L) is calculated numeri- 
cally for different values of the parameter being varied to 
effect the transition and fit to Eqn. [27] The value of the 
parameter which gives the best fit is where the transition 
takes place. In Fig.(j4|) we show that the BKT transi- 
tion point from CSF to CMI is at l/J T = 0.887(1) for 
J± = 1. This method has been used to successfully locate 
the thermal BKT transition for a single layer XY model, 
withal or without frustration 41 -, and has also been used 
to detect non-BKT thermal transitions driven by half- 
vortices (which yield A — 8/ir) in spinor condensates 4 ^. 
Using this technique for several values of J± we obtain 
the boundary for the CSF-CMI transition as shown in 
Fig.H 

The Ising transition requires a different method for its 
detection. Since the chiral order is truly long ranged, we 
can use the method of Binder cumulants to accurately 
locate the transition. For an order parameter to, the 




FIG. 4. Helicity modulus T versus 1/J T for different system 
sizes for J± — 1 as also obtained in ref. fl3l . (Inset) RMS error 
of fit to the BKT finite size scaling form of V shows a deep 
minimum at the transition, at 1/J T =0.887(1), and yields a 
jump AT^ 0.637, close to the BKT value of 2/V. 

Binder cumulant is defined as 

For our system, the order parameter is given by 

m =7iE(- 1 ) lJ -' ( 29 ) 

ir 

where Ji T is the current around a plaquette normal to the 
time direction, i and r are respectively the coordinates in 
the || and time directions. If a(i,r), a(i-\-l,r), b(i + l,r) 
and b{i, r) are the vertices of the plaquette going around 
clockwise, 

J lT = J,, [sin (tpf +l T - <pl r ) + sin (ipl T - <A b +i,r)] 

+ Jj_ [sin (^ +1 , r - ^ +1>T ) + sin « T - ^ T )](30) 

The order parameter given by Eqns. [2FJ and [3D] is like 
a staggered magnetization generated by the loop cur- 
rents. Since the current in a given loop can be clockwise 
or counterclockwise, this order parameter is of the Ising 
type. The values of Bl calculated for different L are 
equal only at the fixed points of the system. In addition 
to the low and high temperature fixed point, they will 
be equal at the location of a continuous transition sepa- 
rating an ordered phase from a disordered phase. Thus, 
plotting Bl as a function of the tuning parameters for 
different values of L reveals the location of the transition 
at the point where the curves intersect as shown in Fig. 
[5] where we plot Bl as a function of 1/J T for J± = 1. 
Curves for different L intersect at 1/J T = 0.981(4) which 
is the transition point from CMI to MI phase. 

It is important to note that the method of Binder cu- 
mulants does not assume any particular universality class 
for the transition. To verify that the transition in ques- 
tion is indeed of the Ising universality class, as expected 



T 




4 1 1 1 1 1 1 

0.96 0.98 1 1.02 

FIG. 5. Binder cumulants for the staggered current versus 
1/J r (for different L for J± — l) intersecting at a continuous 
transition at 1/JV =0.981(4) as also obtained in ref . [ilil. (In- 
set) Critical susceptibility versus L gives the ratio of critical 
exponents 7/f«1.72, very close to 2D Ising value 7/^ = 7/4. 



on the grounds that the CMI breaks a two-fold discrete 
translational symmetry (or time reversal symmetry, we 
have calculated the susceptibility \ for the onset of chi- 
ral order. At the transition, x is expected to diverge 
with system size as x ~ L 1 ^, where 7 and v are re- 
spectively, the conventional susceptibility and correlation 
length critical exponents. For the classical 2D Ising uni- 
versality class, ■y/v = 7/4, while we find 7/1/ « 1.72(inset 
of FigJS]), which demonstrates that our transition is in- 
deed in the Ising universality class. By performing this 
proceedure for several values of J± = 1 we obtain the 
Ising transition line as shown in the Fig. [3l 

In the previous studies of .Hxy it has been shown 
that there exists only one transition when J± = Jm and 
for large anisotropy the system exhibits two separate 
transitionSi 43 i 44 In our study, we find the CMI phase at 
the fully symmetric point Jj_ = J|| = J T indicating two 
transitions. The absence of the second phase transition 
in the previous study could be due to the small system 
sizes considered. 



V. DENSITY MATRIX RENORMALIZATION 
GROUP STUDY 

We have performed a finite size density matrix renor- 
malization group (FS-DMRG) study on the model de- 
scribed by Eqn.[T]to understand different quantum phase 
transitions and compute an accurate phase diagram i 45 i 46 
In our FS-DMRG calculation we have mapped the 2-leg 
Bosc-laddcr into a single chain with appropriate hopping 
elements. Calculations were performed up to a system 
length of 200 sites (which is equivalent to 100 rungs of 
the ladder system) at a filling one boson per site. We 
keep six states per site and retain 200 density matrix 
states in our calculation. The error due to the weight 




FIG. 6. Phase diagram of the FFBH model in Eqn. [T] ob- 
tained using DMRG in ref. [l3l . Both the classical and quan- 
tum models exhibit a chiral Mott insulator (CMI) at inter- 
mediate correlations, intervening between a chiral superfluid 
(CSF) and an ordinary Mott insulator (MI). 

of the discarded states is expected to be less than 10 ~ 5 . 
By calculating different relevant physical quantities we 
obtain an accurate ground state phase diagram of Eqn. Q] 
as explained below. 

We analyze the CSF-CMI and CMI-MI transitions us- 
ing a finite-size analysis of the momentum distribution 
function and the rung current structure factor respec- 
tively. These are defined as 

< k ) = \ E elk(x ~ x ' ] [< a * a ^ + ^] ( 31 ) 

x,x f 

and 

5 #) = JlE elfc( ^'W}, (32) 

where j x = i {a\b x -b\.a x ). 
At a BKT transition, 

(ala x ,) = (blb x ,)~l/\x-x'\ 1 l\ (33) 

so n(0)L~ a with a = 3/4 is independent of the system 
length L. Thus curves of n(0)L -3 / 4 as functions of the 
tuning parameter for different L will intersect at a point 
which we have used to locate the transition as shown in 
Fig. [Vl In the top panel of Fig. \V\ we plot n(k = 0)L~ 3 / 4 
versus U jt with tj_ = t. The curves for different L inter- 
sect at Ud/t m 3.98(1) showing the CSF-CMI transition. 
The inset shows that the charge gap in the insulator also 
becomes nonzero at this point, thus providing a nontriv- 
ial consistency check that we have correctly located the 
superfluid to insulator transition. 

In order to see that such a intersection is not obtained 
for a < 3/4 (corresponding to real space correlations 
decaying with a power law faster than l/\x — x'] 1 / 4 ), we 
plot, in the bottom left panel of Fig. [Vj n(0)L~ a with 



a = 0.6 and find that the crossing point of these curves 
for different system sizes drifts significantly instead of 
coinciding. For example, the crossing point between L = 
70 and L = 90 is at U/t = 4.040, where as it shifts to 
U/t = 4.048 between L = 60 and L = 70. Within the 
CSF state, however, we expect real space correlations to 



decay slower than 1 / 



,/|l/4 



so that curves of n(0)L~ 



with a > 3/4 are simply expected to cross at smaller U 
(i.e., deeper into the CSF and away from the CSF-CMI 
transition). This is also shown in the bottom panels of 
Fig. |V1 where the crossing point is much better defined 
for a = 0.75 (bottom middle panel, critical point) and 
a = 0.8 (bottom right panel, CSF state). 

To summarize, the study of the momentum distribu- 
tion shows that the ground state of the FFBH model has 
critical supcrfiuid correlations, with a > 3/4, but that 
there is no such critical supcrfiuid for a < 3/4, suggesting 
that a = 3/4 is the endpoint of a line of critical points. 
Finally, the CSF-CMI transition located in this manner 
from the momentum distribution completely agrees with 
the superfluid to insulator transition inferred from the 
onset of a charge gap shown in the inset of Fig. [Vj 
Such a careful analysis thus unambiguously shows that 
CSF-CMI transition is of the usual BKT type. 

Having located the CSF-CMI transition, we next turn 
to the CMI-MI transition. Based on the classical model 
study, the CMI-MI transition is expected to be an Ising 
transition. In order to locate this transition for the FFBH 
model using DMRG, we use the scaling form of the rung 
current structure factor Sj(n). Since the chiral order is 
a staggered pattern of current loops, the ordering vector 
is at k = it. To obtain the critical point we use the full 
scaling form 



Si(*)LW v = F (([/ -U c2 )L^ 



(34) 



where /3 and v are respectively the order parameter and 
correlation length critical exponents, and U is the tuning 
parameter with the transition located at U C 2- For the 
classical 2D Ising universality class, (3 = 1/8 and v = 1. 
As a result of the above scaling, curves of Sj(ir)L 2 PI 1 ' 
for different system sizes intersect at the transition point 
U C 2 ~ 4.08(1)2 as shown in Fig. [5] Further, with this 
choice of exponents to rescale the structure factor and 
distance to the transition, we obtain a complete data 
collapse as shown in the inset of Fig. [5J These are un- 
ambiguous indicators that we have correctly located the 
CMI-MI transition. The phase diagram obtained using 
this technique is shown in Fig. [5] 

Finally, in order to show that the CMI is a charge 
gapped state of bosonic matter, we have computed the 
charge gap for t = t± = 1. The thermodynamic extrap- 
olation of this charge gap is shown in the inset of the 
top panel in Fig. [V] The error on the gap is w 0.01. 
Again, this unambiguously points to a gapless CSF state 
for U < U c i, while the states identified as CMI and MI 
have a nonzero charge gap. 

Using DMRG, we have mapped out the momentum 
distribution in the various phases - CSF, CMI, and MI. 




3.90 
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FIG. 7. Top panel: n(0)ZT 3/4 as a function of U/t for different 
system sizes, showing a crossing point at U c .i/t = 3.98(2) 
which we identify as the CSF-CMI transition point which is 
in the BKT universality class. The inset shows the charge 
gap opening up at this transition into the insulator. Bottom 
panel: n(0)L~ a as a function of U/t for three different values 
of a = 0.6 (left), 0.75 (middle) and 0.8 (right) and different 
system sizes. The curves can be seen to cross sharply for 
a = 0.75 and 0.8 whereas they do not for a = 0.6. 



An analytical discussion of this momentum distribution 
in the CSF state is given in Appendix B. In all the phases, 
we find two peaks at k = 0, 7T, as seen from the DMRG 
results shown in FigJU While the two peaks in the CSF 
phase are sharp and grow rapidly with system size, a 
careful scaling analysis shows that the peaks in the mo- 
mentum distribution do not grow with system size in the 
CMI and MI. While this is hard to detect by eye in the 
CMI, which shows a sharp peak even for L = 100, the 
two peaks are clearly extremely broad in the MI. 



VI. PHYSICAL PICTURES FOR THE CMI 

Having identified the remarkable intermediate CMI 
state using careful numerical studies, we next turn to 
simple physical pictures for this state. We first describe 
the CMI as an exciton condensate approaching it from 
the Mott insulator. Next, we argue, from the superfluid 
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FIG. 8. Rung current structure factor Sj(Ti)L 213 ^ 1 ' versus 
U/t at t± = 1. The intersection point yields the CMI-MI 
Ising transition at U C 2 ~ 4.08(l)t as also obtained in ref . Il3l. 
Inset shows 5 , i (7r)L 2/3/ " versus 5L 1/v with 5 = (U - U c2 )/t, 
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FIG. 9. Momentum distribution n(k), defined in Eg. 1311 com- 
puted using DMRG, plotted in the CSF, CMI, and MI regimes 
of the phase diagram showing symmetric peaks at k = and 
k — tv. While the CSF exhibits sharp peaks which grow with 
system size, the CMI and MI phases have finite peaks which 
do not scale with system size. 



side, that the CMI may be alternatively viewed as a vor- 
tex supcrsolid. 



A. Exciton Condensate 

Consider the regular Mott insulator state with an in- 
teger uq > 1 number of bosons at each site at strong 

I 



repulsion U 3> t,t±. The low energy excitations about 
this state correspond to 'doublons', obtained by having 
an extra particle at a particular site leading to an occu- 
pancy no + 1, and 'holes', obtained by removing a particle 
leading to an occupancy no — 1. The energy cost of a cre- 
ating a doublon is Uno—^Jb, while the energy cost to create 
a hole is fi — U{n$ — 1). Once we have a single doublon or 
hole, this can move around by boson hopping - the rel- 
evant hopping matrix element for these particles can be 
shown to be just the bare boson hopping amplitude (±i, 
or t±) enhanced by a factor of (no + 1) for doublons and 
no for holes. The resulting low energy bands of doublons 
and holes thus have energies 

£ d (k) = Un -fi-(n + l)E k , (35) 
Efc(fc) = n - U(n - 1) - n E k , (36) 

where — E k is the lowest band dispersion of a free 
sing le particle on the ladder, given by —Ek = 
— ■\/{2t cos k) 2 + 1\, which has degenerate minima at 
k = 0,7r. We can pick \i such that we have particle- 
hole symmetry in the sense that the energy required to 
create the lowest energy hole (at k = 0, tt) is the same 
as that required to create the lowest energy doublon (at 
k = 0,7r). This leads to fj, = U(n - 1/2) - E /2, and 



^(0/tt) = e h (0/ir) = | - (n + l -)E . (37) 



At large U ^> t, t±, this energy is positive, so that 
doublons and holes are suppressed in the Mott insu- 
lating ground state. A crude estimate of the transi- 
tion from the Mott insulator to a superfluid phase is 
given by demanding that £d,h(fi/Tr) = 0, which yields 
U c (n ) = (2n + l)yjtt 2 + t\, which for t±/t = 1 and 
no = 1 yields U c = 6.7i, within a factor-of-two of the nu- 
merical DMRG result. More importantly, this calculation 
shows that the excitations of the regular Mott insulator 
are similar to gapped particles and holes in a semicon- 
ductor, and the Mott insulator to superfluid transition 
is similar to metallizing a semiconductor by reducing its 
gap. This suggests the possibility of forming doublon- 
hole bound states, akin to excitons, in the vicinity of the 
insulator to superfluid transition where the gap is small. 
Since our "Mott semiconductor" has multiple minima in 
its band dispersion, both 'direct excitons' and 'indirect 
excitons' are possible depending on whether the doublon 
and the hole are at the same momentum (both at k = 
or both at k = 7r), or are at different momenta (one at 
k = and the other at k = ir). 

Insight into the nature of the CMI is obtained by not- 
ing that the rung current operator j a b(%) = —it±(al,b x — 
b\.a x ) can be re-expressed in the Mott regime as 



j ab {x) |Matt)=-itx Vno (no + 1) (4 (*) h\ (x) - d\ (x)hl (x)) |Mott) 



(38) 



where the operators $ a , b (x) and h', b (x) create doublons 
or holes on leg a/6 on the rung labelled by x, and |Mott) 
denotes a caricature of a Mott state with precisely no 
bosons at each site. Focusing on the low energy doublon 
and hole modes amounts to projecting these creation op- 
erators to the lowest dispersing band, and to the vicinity 
of k = 0, 7T. This yields 

h\(x) = u ~hl(x) + v {-l) x ht(x) (39) 

h\{x) = v hl(x) + u (-l) x hUx) (40) 

4(a) = u Q dl(x) + v (-l) x dt(x) (41) 

4(x) = v dl(x) + uo(-l) x di(x) (42) 

where h) , are the low energy creation operators in the 
vicinity of the indicated momenta. In terms of these, the 
rung current operator acting on the Mott state becomes 

j Qb (x)|Mott) = -it ± {ul-vl){-lf{dlhl - d\~hl)\Mott). 

(43) 

This shows that the current operator behaves as a com- 
posite operator that creates an 'indirect exciton'. It 
suggests that the CMI state, in which the current op- 
erator on the rungs has long range staggered order, 
may be viewed as an 'indirect exciton' condensate,— ob- 
tained by condensing this composite operator, to yield, 
{d\h\ — d\h\) ~ i 1 ^, in the absence of any doublon/hole 
condensation, i.e., with {d Q ^) = (h^^) = 0. We note 
that the preceding physical insight into the CMI phase 
does not yield a simple energetic reason for why there 
should be such an intermediate phase between the regu- 
lar Mott insulator and the chiral superfluid. An alterna- 
tive scenario, within Landau theory, is that the combined 
presence of soft hole/doublon modes and such low energy 
excitons might render the transition first order with no 
intervening CMI state. 



B. Vortex Supersolid 

The CSF which possesses staggered loop currents can 
be understood as a vortex-antivortex crystal. In such 
a situation, the vortices and antivortices are generated 
due to the effect of frustration, and due to the repulsion 
between them, they arrange themselves in and antiferro- 
magnetic order. When U is large, the vortices delocalize 
to form a vortex superfluid. The intersting thing hap- 
pens when there exists a small number delocalized and 
coherent defects (vacancies or interstitials) in the vortex 
crystal. These can destroy the superfluidity but while 
preserving the underlying vortex crystal structure. The 
resulting state can be regarded as a vortex supersolid, 
and it is equivalent to the CMI in our model. 
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FIG. 10. Boson correlator G^ 2 \x — x') = (a^a^i) computed 
using variational Monte Carlo method for: (i) the mean field 
CSF state (circles) with off-diagonal long range order, (ii) a 
correlated CSF with a Gutzwiller factor (squares), and (iii) 
the CMI state with a long-range Jastrow (triangles) which 
leads to an exponential decay of (x — x 1 ). The oscillations 
stem from low energy boson modes at momenta k — and 
k — tt. (Inset) Nonzero staggered current order ( — l) x (Ju(x)) 
in the mean field CSF and in the CMI. 



VII. VARIATIONAL WAVEFUNCTION FOR 
THE CMI 

In order to go beyond mean field theory within a wave- 
function approach, one can incorporate Jastrow factors 
which build in interparticlc correlations. Such Jastrow 
wavefunctions for an N-boson system take the form 

*(n, r 2 , . . . ,rjv) = e" ^ ^ ri ~ ri) ^ MF {ru r 2 , ...,r N ) 

(44) 

The simplest example of such a state is the Gutzwiller 
wavefunction for which v(ri — rj) is nonzero and positive 
for ri = Vj and is zero if ri ^ Tj . Such a Gutzwiller factor 
builds in correlation effects by suppressing those config- 
urations in which multiple bosons occupy the same site. 
Such a Gutzwiller wavefunction is however inadequate 
to describe superfluid and insulating states of bosons. 
The reason is that such short range Jastrow factors do 
not correctly encode the small momentum behavior of 
the equal time density structure factor. For a superfluid 
which supports a linearly dispersing sound mode at low 
momentum, we need to have the Fourier transformed Jas- 
trow factor v(q — > 0) ~ 1/g, while a gapped insulator in 
ID should have v(q — > 0) ~ 1/g 2 , as discussed in Ref. 26. 

In order to describe the chiral Mott insulator, we there- 
fore fix ^mf to be the mean field chiral superfluid, and 
choose the Jastrow factor to have two pieces, a short 
range (on-site) Gutzwiller piece vsr and a long range 
part vlr- We parametrize the on-site Gutzwiller factor 
to have strength g$, and choose the long range Jastrow 
factor to have a Fourier transform, VLR(q) = 1/(1— cosg), 
and with a strength g^. Using this paramctrization, we 
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FIG. 11. Behavior of the density structure factor S(q) in 
the chiral Mott state with both short and long range Jastrow 
factors, showing 5(g) ~ q 2 at small momenta. Inset depicts 
the singular long range Jastrow factor vlr ~ 1/q 2 which is 
necessary to obtain this insulating phase. 



have computed the properties of the resulting variational 
wavefunction using standard Monte Carlo sampling, av- 
eraging over 10 6 -10 7 boson configurations on the two-leg 
ladder for system sizes upto L = 500. 

Fig. [10] shows the computed 2-point correlation func- 
tion G^ 2 \x — x') = (a\.a x ,) for different choices of gs 
and gL- For gs = Ql = 0, we find the mean field result 
where (r) saturates to a nonzero value, while exhibit- 
ing oscillations which correspond to the fact that there 
are boson modes at both momenta k = and k = ir 
making up the condensate wavefunction. For gs = 1 and 
gi = 0, the local Gutzwiller factor suppresses multiple 
occupancy - this weakens but does not destroy the off- 
diagonal long range order present in the mean field state. 
For gs = 1 and = 0.1, we find a rapid exponential de- 
cay of G^'ix — x') (shown in Fig. 8), with a correlation 
length £ ~ 25 lattice spacings obtained from a fit to the 
numerical data. This state thus has exponentially decay- 
ing boson correlations indicating a gapped ground state. 

This is in agreement with the behavior of the boson 
density structure factor S(q). For gs — 1 and gi = 0.1, 
we find, as shown in Fig. II H that the structure factor 
exhibits a ~ q 2 behavior at small momentum which sug- 
gests, within the Feynman single-mode approximation, a 
gapped ground state. From the inset of Fig. [TJJ we find 
that this state nevertheless has an average nonzero stag- 
gered current in each plaquette which is inherited from 
the mean field state, although it is slightly suppressed 
due to the on-site Gutzwiller correlations. This Jastrow 
correlated state thus provides a simple variational ansatz 
for the Chiral Mott insulator ground state. 



VIII. EXPERIMENTAL DETECTION 

How might the CMI be realized in experiments and 
detected? One possibility would to try to realise it in a 
Josephson junction ladder with a flux of hc/4e per pla- 
quette. The detection of this phase would then involve a 
measurement of the resistivity to verify that it is an in- 
sulator along with a measurement of local magnetic mo- 
ments to detect the pattern of staggered moments gen- 
erated by the circulating currents. An estimate shows 
that junctions with a coupling of about 1 Kelvin and lat- 
tice parameter of 10 (im could produce staggered fields 
of about 1 nT resulting from the staggered moments, 
which can be measured by high precision Supercondcut- 
ing Quantum Interference Device (SQUID) microscopy. 

Another possibility is to use synthetic magnetic fields 
to generate the tt flux in a system of cold atoms in an 
optical lattice. The n flux will then produce twin peaks 
in the single particle momentum distribution, which will 
be sharp in the CSF and get broadened by the interac- 
tion in the CMI and MI. There are two possible ways to 
distinguish between the CMI and the MI experimentally 
1) Intermode coherence can be tested for by Bragg rein- 
terference of the k = and k = tt peaks4£ The CMI will 
display coherence between the two modes while the MI 
will not. 2) The presence of staggered current patterns 
can be detected directly using a quench experiment.— 
The staggered currents exist on bonds of a lattice cir- 
culating in a particular way (clockwise or anticlockwise) 
around the plaquettes. If the hopping on certain bonds 
(say those between the legs of the ladder) were to be sud- 
denly turned off by a quench, the current loops would be 
broken resulting in charge accumulation and depletion at 
lattice sites at later times, resulting in a time-dependent 
staggered modulation of the local charge density. This 
density modulation can then be used to back out the 
original pattern of currents flowing on the bonds of the 
lattice. Such modulations would also be present in the 
CSF state, where the density modulations and oscilla- 
tions would be likely to persist for longer, whereas they 
would be short lived in the CMI due to heating effects 
from the quench- induced nonequilibrium currents. The 
observation of such quench induced density dynamics (al- 
beit short-lived) together with a charge gap would be an 
unambiguous signature of an intermediate CMI state. 



IX. GENERALIZATIONS 

We next discuss potential generalizations of our work 
along various directions - although we do not have, at 
this stage, detailed numerical studies on these examples, 
they are sufficiently well motivated physically and worth 
examining in more detail in the future. We begin by 
noting the generalization to the case of staggered flux on 
the ladder, where we have alternating fluxes $, — $ with 
< $ < 7r. We then discuss a possible generalization of 
the FFBH model to 2D. 
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FIG. 12. Left panel: Sketch of a staggered loop current state 
in a 2D CMI state in a fully frustrated Bose Hubbard model 
of weakly coupled two-leg ladders with 7r-flux per plaquette. 
Right panel: Sketch of the [11] edge depicting edge current 
along such a 2D CMI state. 



A. Staggered flux on a 2-leg ladder 

Instead of a uniform 7r-flux per plaquette, let us imag- 
ine having staggered fluxes $, — $ on the two-leg ladder, 
with < $ < 7T. The resulting band dispersion of nonin- 
teracting particles exhibits a single non-degenerate min- 
imum at k = 0. Condensing into this minimum leads 
to a superfluid with a unique current pattern, which al- 
ternates from plaquette to plaquette but whose sense is 
entirely locked to the underlying flux pattern since the 
Hamiltonian itself breaks time- reversal symmetry. In this 
sense, moving away from the fully frustrated limit to- 
wards a staggered flux state is like turning on a "magnetic 
field" which couples to the underlying Ising order asso- 
ciated with loop currents. Thus, the staggered loop cur- 
rents survive to arbitrarily large U /t, in much the same 
way as the paramagnetic phase of the usual Ising ferro- 
magnet acquires a nonzero magnetization in a uniform 
magnetic field at any finite temperature. Turning to the 
quantum phase diagram of the staggered flux Bose Hub- 
bard or quantum XY models, we lose the distinction be- 
tween the CMI and the MI, so the Ising transition associ- 
ated with this distinction will be replaced by a crossover. 
We thus expect to have a single BKT transition from a 
superfluid to an insulator, with non-vanishing staggered 
loop currents at all values oiU/t. Numerical studies con- 
firming this prediction would be valuable. 



B. 2D analogue of the bosonic CMI state 

The 2-leg ladder FFBH model we have studied here 
provides a nice example of an unusual weak Mott insula- 
tor state which supports loop currents since the virtual 
boson number fluctuations permit the bosons to 'sense' 
the flux on each plaquette even though they are insu- 
lating at longer length scales. Imagine weakly coupling 
many such ladders lying next to each other to form an 
anisotropic version of the 2D FFBH model, where each 
square plaquette still has 7r-flux, but where the boson 
hopping takes on intra-chain values ±t along alternate 



chains, with a modulated interchain hopping which takes 
on a value t± within a ladder (intra-ladder rung hopping) 
and t'j_ between adjacent ladders (inter-ladder hopping). 
The sketch of this model is shown in Fig. [T2] (left panel) . 
When t' ± = 0, we get a set of decoupled FFBH ladders 
which support a fully gapped CMI state. We thus ex- 
pect that this phase will be stable so long as the gap 
in the CMI state is much larger than the inter-ladder 
hopping The charge gap in the CMI would render 
this stable against Bose condensation into a 2D super- 
fluid. Similarly, the Ising gap associated with the broken 
Z% time-reversal symmetry would render this 2D CMI 
stable against disordering of the current loop order for 
small values of t' ± . Based on the fact that vortices repel 
each other, and would prefer a global (-7T, 7t) antiferro- 
magnetic ordering pattern in 2D, we expect the current 
loops on the different ladders to phase lock in-step pro- 
ducing a fully gapped 2D CMI state with staggered loop 
currents as shown in the left panel of Fig. [T2] Such a 
CMI state would lead to a chiral edge current along any 
of the [11] edges of the lattice as shown in the right panel 
of Fig. [T2J Numerical studies of weakly coupled ladders 
within a classical XY model calculation would be valu- 
able to explore the full phase diagram as a function of 
U/t, tj_/t, and t'Jt. 



X. SUMMARY 

We have computed the accurate phase diagram of the 
FFBH model using two different numerical techniques: 
Monte Carlo simulations and the FS-DMRG method. We 
have consistently obtained qualitatively similar phase di- 
agrams depicting three different quantum phases: CSF, 
CMI and MI. The CMI phase is sandwiched between 
the CSF and the MI phase, and is a remarkable state 
of bosonic matter, being an insulator but with staggered 
loop currents which spontaneously break time-reversal 
symmetry. We have confirmed that the transition from 
CSF-CMI is BKT type and the transition from CMI to 
MI is of Ising type. We have also shown that the boson 
momentum distribution has different properties in these 
phases. We present the two different physical pictures for 
the CMI as an exciton condensate and the vortex super- 
solid. A variational wavefunction which can explain the 
CMI phase is discussed. Finally, we have proposed pos- 
sible experimental signature of different quantum phases 
which can be realizable in the experiments involving cold 
atoms in optical lattices and Josephson junction arrays, 
and discussed higher dimensional realizations of the CMI 
in a FFBH model. 
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Appendix A: Derivation of classical XY model 

We define the new rotor variables to satisfy the follow- 
ing commutators: 







ct,c p x , =o 



(Al) 
(A2) 
(A3) 
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where a/j3 = a,b. In terms of these, the Hamiltonian 
takes the form 

F rotor = -2tJ2 cos(v>£ - vl+i) + 2*5Z cos ^ - 

a; x 

-2tj.5^cos(^-^) 

+ 1 E [( £ S) 2 + (4) 2 )] - mE( £ " + 4), (A4) 

where we have allowed the rotor hopping to be propor- 
tional to the original boson hopping, with tj_/t = tj_/t. 
We expect the proportionality constant to be such that 
t = tip 2 ~ tp where p is the boson density per site. 

At p — 0, which corresponds to an average of one 
boson per site, we can easily go to the path integral 
representation. We start with the partition function 
Z = Tr(e~P Hrotat ) and use Trotter discretization to set 
cxp(-^iJ rot or) = [exp(-eif rotor )] iT , with e = f3/L T . The 
partition function then takes the form 

Z= E ({%,o}|e- eiW \W x ,l t -x})... 
{W}} 

...({<^x,i}|e- eH — K^.o}) (A5) 

where we have chosen to compute the trace in the angle 
basis and intermediate states are also in this basis. Here, 
single curly brackets refer to an angle-configuration at 
a fixed 'time'-slice, while double curly brackets denote 
angle-configurations over all of space- 'time'. Let us con- 
sider the matrix element 

({y x , T }\e~ eH * ot °*\{v x ,T-i})- (A6) 



with 2et = J||, 2et± = J±, and 1/eU = J T - We see 
that this has the form of an anisotropic XY model. In 
order to get the properties of the quantum model at a 
fixed inverse temperature fit, we must take the 'time'- 
continuum limit of e —> 0, sending Jh — s- 0, J± —> 0, and 
J r — > oo, while keeping fixed J±/J\\ =t±/t and J T Jn = 
2t/U. The inverse temperature fit is then given by etL T 
and thus depends on the chosen value of et (which must 
be taken to be very small) and the size of the simulation 
cell in the 'time'-direction. We set e = \j\j2Ut which 

leads to J\\ = J T = y/2t/U. 

Appendix B: Momentum Distribution 

The main result of this appendix is to show that the 
peaks in the CSF phase of the ladder model have a non- 



Sincc e — > 0, we can set 

(W}|e- £iW |{^, T -i}) 

« ({^, T }|e- eV —e- eK — \{<p x , T -i}) 
« e- eXW «^-»({%,r}|e- eif — |{Vx,r-i}) (A7) 

where V and K are the potential ((^-diagonal) and kinetic 
((^-off-diagonal) terms in the rotor Hamiltonian. Explic- 

itly, 

Kotor(i>x,r}) = -2t5^COs(^ lT ~ VS+l.r) 

X 

+ 2*5^oos(^ iT -^ +liT ) 

X 

-2i ± ^cos« T -< T ). (A8) 

X 

and 

i^rotor = -^^[cOs(^ r -^^ T+1 )+COs(^. r -^ jT+1 )] 

(A9) 

This gives us the partition function in terms of a classical 
action in one higher dimension 

Z = E e- s c\ +1 M (A10) 

where 



(All) 



I 

universal power law divergences, with equal exponents, 
at k = 0, 7r. This is consistent with our DMRG simu- 
lation results of the FFBH model which show two such 
equal height and power law divergent peaks in momen- 
tum distribution in the CSF state (shown in Fig. |H] and 
discussed in Sec. V). 



To analytically calculate the momentum distribution 
in the CSF on the ladder, we begin by carrying out a 
classical minimization of the rotor potential energy in 
Eq. IA4| which leads to the equations 



I 

Scl +1 = - E [ J H C0S (^+l,r - <r) - h COB(^ +liT - < T ) + J± COs(^ |T - < T )] 

XT 

- J T [o»(^ iT+1 ~ <p% tT ) + C0S(^. T+1 - <p b XiT )\ , 

XT 
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t[sin(^ - v a x+ 
t[$m{(p b x - tp b x+1 ) +sm(cp\ 



,6 



■<Pl 



It is straightforward to check that we satisfy these equa- 
tions by substituting cos</5° = u , cos(p x = vq, sin^" = 
±uo(— l) x , and sm(p x = ±ito(— l) x , where uo,vq are de- 
fined via Eqns l4l5l For now, let us stick with the '+' 
solution. We can write this solution as tp x = (— l) x ty, 
ip\ = (— l) x (n/2 — *) where cos* = w and sin* = v . 
To take small fluctuations into account, we must set 



= (-1)^/2 + 



Vx 

J J 



(B3) 
(B4) 



rrsmall 2 

-"rotor ~~ C 



cos(2vI/)^[(^-^ +1 ) 2 



' f ± sin(^-^) 



</>*)= 0, 

: 0. 



(Bl) 
(B2) 



U 



r 



Using the above minimization condition, we find the 
Hamiltonian for small fluctuations, 



(B5) 



Further, for small fluctuations, we must set exp(— icj)) « 
l—i(f>, so that the full commutation relations are replaced 
by 



i5 afi 6 XiX , 







(B6) 
(B7) 
(B8) 



so that and C have exactly the same commutation rela- 
tions as position/momentum variables. Let us now define 

I 



I 

a shifted angular momentum C" = — fi/U to absorb 
the chemical potential term, which leads to the Hamilto- 
nian (upto an additive constant we will ignore), 



rrsmall j. 



r cos(2*)^[(C-C+i) 

X 

+ Usin(2*)^(C-^) 2 

X 

x;[ocs) a +(4) a )' 

a; 

Fourier transforming, we find 



U 
~2 



flJSto=5 E [2tcos(2*)(2-2cos«) + 2t X 8in(2*)]^ 
- 1 £ 2fa sm(2*)(^#_,+0° + £ ^ £°£° 9 



(B9) 



(BIO) 



This is like a problem of two coupled oscillators for each 
q, with each oscillator having a 'mass' = 1/U, but with 
a 2 x 2 spring constant matrix having eigenvalues 

K+{q) = 2tcos(2*)(2-2cosg) (Bll) 
K-(q) = 2icos(2*)(2-2cosg) + 4t. ± sin(2*).(B12) 

This leads to a collective mode spectrum with two 
branches, uj± = \/UK±{q), where uj+(q) = cq for small 
q behaves like an 'acoustic phonon' mode, while uj~(q) is 



I 

like a gapped 'optical phonon' mode. Using the fact that 
sin(2*) = i ± /^P ± +it 2 and cos(2*) = 2t/Jp[+4P, 
we find that the 'speed of sound' is 



AUt 2 



AUpt 2 



(B13) 



while the minimum gap to the 'optical phonon' mode is 
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at q = and given by 
w-(« = 0) = 



AUpt] 



(B14) 



These should be viewed as analogs of a Bogoliubov theory 
result for a ID system with no Bose condensate (i.e., no 
ODLRO). The off-diagonal one-body density matrix may 
be computed by noting that 



(<4«y) 



»*A_(x, 1 /)/i(^*-^;)\ 



(B15) 



where A±(x,y) = (— l) x ± (— 1) J/ . Within the Gaussian 
theory, we get 



(B16) 



The expectation value in the exponential can be evalu- 
ated as 

W x ~ 4>y) 2 ) = W%) [2 - 2cos(g(x - y))] 

(B17) 

We find 



P q F- q ) 



i> q 4>- q ) 





/ u 




K-(q) 




1 U 


K + {q) V 


K-{q) 



(pi8) 
(bl9) 



Since K + (q — > 0) ~ q 2 , while K-(q — > 0) remains 
nonzero, the long wavelength (small q) limit of all such 
correlations are given by 



U 



1 T—ql q ^o 



1 



4 V K+(q) 



(B20) 



Using this, and defining 



F{x - y) 



\x - y\ 



(B21) 



we find the long distance behavior of the one-body off- 
diagonal density matrix elements are given by 



(4s 



3 i*A_(x,y) ■ 



pe~ — ^'»'F{x-y) (B22) 
(btb y ) ~ pe-w*-(*.v)+ii A_ (*,v) F ( x _ y ) (B23) 
(4&„) ~ p ^+(*,v)-i%(-l) y F ( x _ y ) ( B24 ) 



Let us expand 



pe 



-i*A+(a;,j/)+if (-1)' 



F(x - y) (B25) 



using which, 



p(cos 2 *+sin 2 ^'(-l)*"*' 

iA-(x, y) sin 2*)F(a; - y) (B26) 



,(fc) = i^(at % )e ife (^) ~ pcos 2 *^fWe ikr 

/ r 

*^^(r)e- Mr e lfcr , (B27) 



L 

+ psin 



2 



where L is the linear system size. We find n aa (k — > 0) 
(pcos 2 <if)\k\~ a and n aa (k -> 7r) ~ (psin 2 *)|& - 7r|~' 
where 



1 

Q = 1 



u 



,1-A>V / ^ (B28) 
4vr V 2tcos2^ 4?r V 4pt 2 



Similarly, n^ik — > 0) ~ (psin 2 \l')|fc| _a and n(,b(fc — > 
7r) ~ (pcos 2 \I>)|fc — 7r| _Q , so that n aa {k) + ribb)(k) ~ 
/9(|fc| _Q + |fc — 7r| 04 ) . We can also calculate cross 
correlators such as n &(fc) = (a\.b k ) and ni, a (k) = 
(&j.a fc ). We find that n a &(fc) = n& a (fc), and that 
«af>(fc ->• 0) ~ (psinvfi cos\l/)|fc|~ Q , while n a ;,(fc -> n) — 
(psin^cos^lfc - 7r|" a . 

We have thus shown that the peaks in the CSF phase 
of the ladder model have a non-universal power law diver- 
gences, with equal exponents, at k = 0, tt. Similarly, two 
(broad) equivalent peaks in the momentum distribution 
also appear in the MI, but they need a strong coupling 
expansion to compute as has been done, together with a 
collective mode analysis, for the 2D FFBH modelii. At 
the BKT transition from the CSF to CMI, an analysis of 
which requires going beyond our harmonic theory, we ex- 
pect the exponent to take on a universal value a = 3/4. 
Such a bosonization analysis will be described elsewhere. 



